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We calculate the two-loop Bethe logarithm correction to atomic energy levels in hydrogen-like 
systems. The two-loop Bethe logarithm is a low-energy quantum electrodynamic (QED) effect 
involving multiple summations over virtual excited atomic states. Although much smaller in absolute 
magnitude than the well-known one-loop Bethe logarithm, the two-loop analog is quite significant 
when compared to the current experimental accuracy of the 1S-2S transition: it contributes -8.19 
and -0.84 kHz for the 15* and the 25 state, respectively. The two-loop Bethe logarithm has been 
the largest unknown correction to the hydrogen Lamb shift to date. Together with the ongoing 
measurement of the proton charge radius at the Paul Scherrer Institute its calculation will bring 
theoretical and experimental accuracy for the Lamb shift in atomic hydrogen to the level of 10 -7 . 
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In 1947 Hans Bethe explained the splitting of 2S , 1 / 2 and 
2Pi/ 2 levels in hydrogen by the presence of the elec- 
tron self- interaction |lj , and expressed it in terms of the 
"Bethe" logarithm. For S states this quantity may be 
represented as a matrix element involving the logarithm 
of the nonrelativistic Hamiltonian of the hydrogen atom. 
In natural units with h = c = eo = 1 and m denoting the 
electron mass, it reads 
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This Bethe logarithm is due to the emission and sub- 
sequent absorption of a single soft virtual photon (it is 
independent of the nuclear charge number Z and depends 
only on the principal quantum number n and the orbital 
angular momentum which is zero for S states) . Over the 
years, QED theory has been developed and refined 0, 
and various additional radiative, relativistic, and com- 
bined corrections have been obtained to face the increas- 
ing precision of the measurements of the hydrogen spec- 
trum 0-Q • These include higher-order relativistic one-, 
two-, and three-loop corrections, nuclear recoil, finite- 
size corrections, and even the nuclear polarizability. The 
modern all-order calculation of the leading one-loop self- 
energy was developed by Mohr in p| and significantly 
improved recently using convergence acceleration tech- 
niques which led to a highly accurate evaluation of the 
fully relativistic Green function Q . One of the conceptu- 
ally most difficult and as well as interesting corrections 
involve nuclear recoil effects. The finite nuclear mass, 
although large as compared to the electron mass, pro- 
hibits the use of the one-body Dirac equation, and alter- 
native approaches such as the Bethe-Salpeter equation or 
nonrelativistic QED have been introduced. Although 
these methods are quite general, no compact formulas 
have been derived for relativistic recoil effects. A few 
years ago, Shabaev tackled the problem of recoil correc- 
tions to hydrogenic energy levels of first order in the mass 



ratio, deriving expressions which are nonperturbative in 
the nuclear charge (see a recent review in and this 
has led to the current highly accurate calculations of rel- 
ativistic recoil corrections. 

Another class of effects, namely binding two-loop cor- 
rections, are quite difficult from a numerical point of 
view. A detailed investigation of these effects has been 
performed only in the last years. The nonperturbative 
treatment (no Za-expansion) of the two-loop bound- 
state corrections has been pursued by various groups in 
Refs. 0, 0, fill H^ . However, these calculations were 
mostly performed for high-Z hydrogen-like atoms. As 
yet, complete results have not been obtained for Z = 1 
(see for example the most recent work in ^3)- ^ n the 
perturbative treatment of the bound-state two-loop self- 
energy correction, one calculates terms in a semi-analytic 
expansion in Za and ln[(Zcv) -2 ]. For S states, the first 
nonvanishing terms read 
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It was a perhaps surprising result, found only a few years 
ago [l^ |. that this expansion has a very slow conver- 
gence: because of the large absolute magnitude of higher- 
order coefficients, many terms have to be included for 
reliable theoretical predictions which in order to match 
the current experimental precision. One of the remain- 
ing unknown but relevant contributions is the two-loop 
Bethe logarithm, which forms the dominant part of the 
problematic nonlogarithmic coefficient Bqq (note that 
theoretical effort in evaluating the one-loop analog of 
this coefficient, Ap,n . has extended over more than three 
decades [3 EH EHIl3|) The two- loop Bethe logarithm 
originates from the emission and absorption of two vir- 
tual soft photons. Using nonrelativistic QED one derives 
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the following expression for this two- loop correction . 
For convenience we set Z = 1, pull out a common pref- 
actor and express remaining in terms of dimensionlcss 
quantities: 
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Here, p k — —id k , G(u>) = — (ff + £*;)] is the nonrela- 
tivistic Green function, E = — l/(2n 2 ) is the Schrodinger 
energy of the reference state with H = p 2 /2 — 1/r, and 
G'(0) = 1/(E — H)' is the reduced Green function with 
the reference state excluded. The w-integrals in Eq. Q 
depend on e\ and e 2 . We are free to choose the relation 
between the parameters e\ and e 2 . Following [l^, we per- 
form the expansion in large e 2 first, and next in large t\. 
All the terms involving \jt\ or l/e 2 are neglected. The 
linear and logarithmic terms, using the same relations, 
have already been considered in 18]. The dependence 
on the e-parameters cancels when the contributions from 
the high- and the low-energy photons are added. The 
constant term, which by definition we call the two-loop 
Bethe logarithm, is calculated numerically here. 

There are two ways to calculate the integrand in 
Eq. (JSJ . The first one relies on the use of known analytic 
expressions for the Schrodingcr-Coulomb Green function, 
which involve the product of Whittaker functions. The 
precise calculation of the w-integrals requires the use of 
large u>\ and u> 2 . This leads to a number of problems, 
including a numerical overflow in the calculation of these 
functions which persists even in quadruple precision, and 
this approach has therefore not been pursued here. In the 
second way, which is chosen here, the Schrodinger Hamil- 
tonian is represented on a numerical grid, as a large sym- 
metric band matrix |l9). Each inversion in Eq. JSJ cor- 
responds to a solution of a linear equation with a known 
right-hand side. This process is quite fast since it scales 
linearly with the number of grid points and is numeri- 
cally stable. For the final evaluation we used 10 5 grid 
points, and we have checked the numerical accuracy of 
the results against those obtained with 2 • 10 5 grid points. 



Having calculated the matrix elements, we proceed to 
the evaluation of the w\- and w 2 -integrals. In accordance 
with the definition of two-loop Bethe logarithm, we first 
fix u}\, integrate over w 2 and expand in large e 2 : 
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As in the case of the one-loop Bethe logarithm, the g 
function finds a representation that is suited for a nu- 
merical computation, 
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with arbitrary M. The a and b coefficients are the first 
terms of the expansion of /(oji,u; 2 ) for large o> 2 at fixed 
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The first coefficient is 
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and the second reads 
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where by by we denote the first-order correction to the 
specified matrix elements. Namely, <f>, E, and H receive 
corrections according to 
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The higher-order coefficients c, d, . . . in Eq. © are ob- 
tained from the fit to the numerical data and are sub- 
sequently used for the analytic integration at large w 2 . 
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TABLE I: Sample values of the g function, defined in Eq. 
for the IS and 2S states. 





Q l S" 


8 (79 c? 





0.000 00 


0.000 00 


5 


-10.28160 


-10.367 94 


20 


-16.560 34 


-16.415 97 


80 


-22.714 02 


-22.439 66 


180 


-26.232 35 


-25.923 09 


320 


-28.699 64 


-28.376 26 


500 


-30.599 22 


-30.268 75 


720 


-32.142 95 


-31.808 43 



The numerical integration over ui2 is performed with a 
well adapted set of 400 grid points, and the accuracy is 
checked by comparison with a calculation involving 200 
grid points. Results of this integration for few chosen 
values of u>i is shown in Table [fl The next step is the 
numerical integration over u>i, 
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In order to perform this integration, one needs to know 
the large-wi asymptotics of g, which for an nS state 
reads, 
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Here, N denotes a nonlogarithmic in a correction to the 
Bethe logarithm induced by a Dirac 5. It has been cal- 
culated in [2H Eq. (12)], and the results for the IS and 
25 states read 



N(1S) = 17.855 672(1), 
N(2S) = 12.032 209(1). 



(15a) 
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The terms proportional to the A, B, C-coefficients in 
Eq. (| X4p . and the omitted higher-order terms are ob- 
tained from the fit to the calculated data. The numerical 
stability of the parameters obtained from the fit, in vari- 
ous ranges of uii , indicates consistency of the numerically 
determined values for g(oJ\) with the analytically derived 
logarithmic terms in Eq. (|14|) . This constitutes a check 



for the large value for the coefficient B$i derived in 0] 
on the basis of the logarithmic asymptotics. 

The integral in Eq. 1)1 3|l is performed in analogy to 
the algorithm presented in Eq. by choosing an ar- 
bitrary value for the parameter M (e.g., M = 1), and 
dropping all linear and logarithmic terms in e\. For oj\ 
larger than 720, the extrapolated values from the fit are 
employed. However, this part of the integral depends 
significantly on the unknown analytic behavior of g(tij) 
at large co, and this is the main source of the integra- 
tion uncertainty. The overall result of the numerical in- 
tegration leads to the following nonlogarithmic terms of 
the order of (a/7r) 2 (Za) 6 /n 3 relative to the electron rest 
mass, whose numerical coefficients we choose to denote 
by b L , 



b L (lS) = -81.4(3) 
b L (2S) = -66.6(3) 



(16a) 
(16b) 



These terms are much larger than the corresponding one- 
loop Bethe logarithms for typical hydrogenic states, but 
suppressed in absolute magnitude by an additional fac- 
tor a (Za) 2 /it. For hydrogen (Z = 1), the above results 
contribute -8.19 and -0.84 kHz to the IS" and 2S states, 
respectively. The other contributions to Bqq are consid- 
ered below; the notation is consistent with that of Ref. 
flij ]. The coefficient Bqq can be represented as the sum 
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The two-loop Bethe logarithm 6l comes from the region 
where both photon momenta are small and has been the 
subject of this work. 6m stems from an integration region 
where one momentum is large ~ to, and the second mo- 
mentum is small. This contribution is given by a Dirac 
S correction to the Bethe logarithm. It has already been 
derived in but not included in the theoretical predic- 
tions for the Lamb shift: 
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bp and bfj originate from a region where both photon 
momenta are large ~ m, and the electron momentum is 
small and large respectively. Finally, &vp is a contribu- 
tion from diagrams that involve a closed fermion loop. 
None of these effects have been calculated as yet. On 
the basis of our experience with the one- and two-loop 
calculations [2lJ we estimate the magnitude of these un- 
calculated terms to be of the order of 15%. This leads to 
the following overall result for the B§o coefficients: 

B 60 {1S) = -61.6(3) ±15%, (19a) 
B 60 {2S) = -53.2(3) ±15%, (19b) 

and to the following corrections to transition frequencies 



Sv(lS) 
5v{2S) 



-6.20(93) kHz, 
-0.67(10) kHz. 



(20a) 
(20b) 
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In the foreseeable future, we may expect to have results 
from a direct numerical calculation of the two-loop self- 
energy at low Z . In addition to improving our knowledge 
of the Lamb shift at low Z, the obtained result will then 
serve as a consistency check between two different ap- 
proaches to bound-state QED. 

With the results obtained in Eqs. I|19all9b[) we are 
in a position to present theoretical predictions for the 
Lamb shifts of IS* and 25 states. Based on the former 
result obtained in and the corrections calculated in 
this work, we obtain: 

v L (\S) = 8172 811(32)(2)kHz, (21a) 
u L (2S) = 1045 005(4) kHz, (21b) 

where the first error comes from the current uncertainty 
in the proton charge radius r p , and the second one is a 
rough estimate of uncalculated terms: bu,byp, as well 
as higher-order two-loop corrections denoted by dots in 
Eq. J2J, and the three- loop binding correction C50 (for a 
recent evaluation of C40 see Ref. [22|). Theoretical pre- 
dictions for the IS state agree well with the experimental 
value of the combined result from the Garching and the 
Paris groups [1 H EsT] : 

f£(lS% xp = 8 172 840(22) kHz. (22) 

Since the uncertainty coming from the proton structure 
dominates theoretical uncertainties for the Lamb shift in 
hydrogen, an experiment at the Paul Scherrer Institute 
(PSI) is currently being pursued to make a precise mea- 
surement of r p from the Lamb shift in muonic hydrogen 
|24j |. This system is much more sensitive to the proton 
charge radius. Once it is measured, the combined hydro- 
gen and muonic Lamb shift will test QED at a precision 
level of 10~ 7 . 
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